Celem poniższych obliczeń jest zbadanie statystyki orbit małych obiektów Ukłądu Słonecznego, znajdujących się głównie pomiędzy orbitami Marsa i Jowisza, znanych jako ,,Pas asteroid’’.

SetDirectory[NotebookDirectory[]];

ast=MinorPlanetData[];

p=MinorPlanetData["Properties"];

ast[[1]]

MinorPlanetData[ast[[1]],"Density"]

periods=N[QuantityMagnitude@MinorPlanetData[ast,"OrbitPeriod"]];

Length[periods]

(*Export["MinorPlanetPeriods.wl",periods]*) (* Niepotrzebne, powyższa linijka załatwia problem powolnego pobierania danych *)

Wygodną jednostką jest okres Jowisza, dominującej, największej (t.j. najmasywniejszej) planety Ukł. Słonecznego, która obejmuje w przybliżeniu ok. 0.001 (Subscript[M, ⊙]), oraz 2/3 masy pozostałych planet.

jp=PlanetData["Jupiter","OrbitPeriod"]//QuantityMagnitude

sp=PlanetData["Saturn","OrbitPeriod"]//QuantityMagnitude

Warto zauważyć, że stosunek okresów orbit Jowisza do Saturna to około 2/5=0.4. Nie jest to przypadek!

jp/sp

periods[[1;;2]]

Na pierwszy rzut oka histogram nie jest ciekawy....

Histogram[periods,{0.2}]

... ale po narysowaniu go z mniejszymi binami, pojawia się charakterystyczna struktura. Położenie przerw i szczytów wydaje się na pierwszy rzut oka przypadkowe.

Histogram[Select[periods[[All]]/jp,0<#<1&],{0.002}]

W tym momencie potrzebne jest trochę wiedzy z analizy matematycznej i matematyki dyskretnej na temat liczb wymiernych. Potrzebujemy ich posortowaną listę, od ,,najmniejszych’’ do ,,największych’’, ale nie w sensie wartości liczbowej, tylko stopnia ,,złożoności’’. Nie ma tu oczywistej kolejności, pomimo że intuicyjnie im większy licznik i mianownik, tym ułamek jest bardziej złożony. Znamy kilka sposobów generowania ułamków, najciekawsze są te, które pozwalają zrobić to bez powtórzeń:

rational=Table[Complement[FareySequence[k+1],FareySequence[k]],{k,1,8}]//Flatten

rational=NestList[1/(Floor[#]+1-FractionalPart[#])&,1,12] (* rational number enumeration *)

rational=Select[NestList[1/(Floor[#]+1-FractionalPart[#])&,1,32] ,#<=1&](* rational number enumeration *)

Ułamki reprezentujące okres planetoid wyrażony w jednostkach okresu Jowisza nakłądamy na histogram.

kirkwood=Histogram[Select[periods[[All]]/jp,0<#<1&],{0.002},
FrameTicks->{rational,Automatic},
PlotRange->{{0,1},Automatic},AspectRatio->1/2,
FrameLabel->{"Okres orbitalny [Okres orbitalny Jowisza]","Liczba asteriod"},Frame->True,
GridLines->{rational,None}]

Po nałożeniu ułamków, obraz wyklarował się. Przerwy i piki są powiązane z ,,małymi’’ w wyżej opisanym sensie ułamkami. Pojawiające się przerwy nazywamy przerwami Kirkwooda. Są one jednym z najbardziej spektakularnych przejawów struktury rezonansowej Układu Słonecznego. Warto zaznaczyć, że potrzebujemy histogramu aby je zobaczyć. Wizualnie, najbardziej  spektakularnym pokazem rezonansów są przerwy w pierścieniu Saturna.

kirkwoodZoom=Histogram[Select[periods[[All]]/jp,1/4<=#<1/2&],{0.001},
FrameTicks->{rational,Automatic},
PlotRange->{{1/4,1/2},Automatic},AspectRatio->1/2,
FrameLabel->{"Okres orbitalny [Okres orbitalny Jowisza]","Liczba asteriod"},Frame->True,
GridLines->{rational,None}]

rational=Table[Complement[FareySequence[k+1],FareySequence[k]],{k,1,8}]//Flatten

gridLinesByIndex=MapIndexed[{#1,(*x‑position*)
Directive[(*style*)
Blue,(*colour*)
Opacity[(1-#2[[1]]/(1+Length[rational]))],(*fade*)
AbsoluteThickness[(1-#2[[1]]/(1+Length[rational]))*2]                  (*thickness*)
]}&,rational];

kirkwood=Histogram[Select[periods[[All]]/jp,0<#<1&],{0.002},
FrameTicks->{rational,Automatic},
PlotRange->{{0,1},Automatic},AspectRatio->1/2,
FrameLabel->{"Okres orbitalny [Okres orbitalny Jowisza]","Liczba asteriod"},Frame->True,
GridLines->{gridLinesByIndex,None}]

Powyższa demonstracja polega na danych z bazy Wolfram-a, które z natury są niekompletne i nieaktualne. Poniżej to samo, ale z użyciem ,,profesjonalnych’’ źródeł danych.

Można spróbować pobrać dane z MPC (Minor Planet Center), trwa to kilka minut.

{czas,raw}=Import["https://minorplanetcenter.net/iau/MPCORB/MPCORB.DAT.gz","GZIP"]//AbsoluteTiming; (* Patience, 188 seconds to download *)

czas

(*raw=Import["data/MPCORB.DAT.gz"];*)

raw[[1]]

raw[[2]]

raw[[3]]

raw[[41;;55]]

head=raw[[41]]

head[[11]]

dat=raw[[44;;]];

dat[[-1]]

a=Table[dat[[ii,11]],{ii,1,10^5}];

Length[dat]

Histogram[a,{0.01}]

Histogram[Select[a^(3/2)/jp,0<#<1&],{0.0005},FrameTicks->{{Automatic,None},{rational,None}},PlotRange->{{0,1/2},Automatic},AspectRatio->1/2,FrameLabel->{"Okres orbitalny [Okres orbitalny Jowisza]","Liczba asteriod"},Frame->True,BaseStyle->FontSize->14,GridLines->{rational,None}]

{czas, nasa}=Import["https://ssd.jpl.nasa.gov/dat/ELEMENTS.NUMBR","Lines"]//AbsoluteTiming; (* Cierpliwości... ... pobieranie trwa 225 sekund *)

czas

nasa[[1]]

nasa[[2]]

nasa[[3]]

nasa[[-1]]

StringTake[nasa[[-1]],{33,41}]

Length[nasa]

aNASA=Table[StringTake[nasa[[ii]],{33,41}]//ToExpression,{ii,3,Length[nasa],1}];

Histogram[Select[aNASA^(3/2)/jp,0<#<1&],{0.0005},FrameTicks->{{Automatic,None},{rational,None}},PlotRange->{{0,1/2},Automatic},AspectRatio->1/2,FrameLabel->{"Okres orbitalny [Okres orbitalny Jowisza]","Liczba asteriod"},Frame->True,BaseStyle->FontSize->14,GridLines->{rational,None},ImageSize->450]

nasa[[3]]

eNASA=Table[StringTake[nasa[[ii]],{42,54}]//ToExpression,{ii,3,Length[nasa],1}];

Histogram[eNASA,{0.005}] (* Mimośrody orbity eliptycznej *)

iNASA=Table[StringTake[nasa[[ii]],{55,62}]//ToExpression,{ii,3,Length[nasa],10}];

Histogram[iNASA,{0.1}] (* Nachylenie płaszczyzny orbity *)